1 Load packages

Load packages needed for data manipulation and analyses.

library(knitr)      #for markdown knitting
library(tidyverse)  #for data wrangling etc
library(colorspace) #for working with colours in visualisations
library(cmdstanr)   #for cmdstan
library(brms)       #for fitting models in STAN
library(rstan)      #for interfacing with STAN
library(ggeffects)  #for partial effects plots
library(DHARMa)     #for residual diagnostics
library(emmeans)    #for marginal means etc
library(tidybayes)  #for more tidying outputs
library(ggridges)   #for creating ridge plots
library(vegan)      #for analysing ecological data
library(EcolUtils)  #for analysing ecological data
library(patchwork)  #for multi-panel figures
library(HDInterval) #for HPD intervals
library(report)

2 Functions

Custom ggplot theme for visualisation.

my.theme <- function(){
  theme_classic() +
    theme(text = element_text(family = "Avenir Next"),
          axis.title.y = element_text(margin = margin(t = 0,r = 20,b = 0,l = 0)),
          axis.title.x = element_text(margin = margin(t = 20,r = 0,b = 0,l = 0)), 
          plot.margin = unit(c(5, 10, 5, 10), units = "mm"),
          strip.background = element_rect(fill = "#CCCCFF"),
          strip.text.x = element_text(size = 20),
          axis.title = element_text(size = 20),
          axis.text = element_text(size = 18),
          legend.text = element_text(size = 15),
          legend.title = element_text(size = 15))
}

3 Load data

reptile_data <- read.csv("reptile_data.csv")

4 Species Richness

4.1 Prepare data

Created an object with coordinates for each plot to assess the effect of latitude.

site <- c(rep("Duval",4), rep("Mourachan",4), rep("Tarcutta",4), rep("Wambiana",4), rep("Undara",4), rep("Rinyirru",4))
site.plot <- c("Duval.DryA","Duval.DryB","Duval.WetA", "Duval.WetB", "Mourachan.DryA", "Mourachan.DryB", "Mourachan.WetA", "Mourachan.WetB", "Tarcutta.DryA", "Tarcutta.DryB", "Tarcutta.WetA", "Tarcutta.WetB", "Wambiana.DryA", "Wambiana.DryB", "Wambiana.WetA", "Wambiana.WetB",  "Undara.DryA", "Undara.DryB", "Undara.WetA", "Undara.WetB", "Rinyirru.DryA", "Rinyirru.DryB", "Rinyirru.WetA", "Rinyirru.WetB")
lat <- c(-30.41734901, -30.40221297, -30.41803398, -30.40110799, -27.77898598, -27.77996097, -27.78425401, -27.77876303, -35.36852497, -35.37876002, -35.36012203, -35.36613797, -20.53070999, -20.526226, -20.53458997, -20.53051997, -18.18705803, -18.24565097, -18.18492902, -18.26552999, -15.05447703, -15.04400703, -15.04891698, -15.04438899)
long <- c(151.622667, 151.624706, 151.61276, 151.629483, 149.032006, 148.981, 149.021332, 148.970948, 147.696893, 147.703341, 147.697579, 147.707788, 146.113426, 146.110754, 146.103876, 146.101471, 144.539753, 144.553799, 144.5324, 144.556613, 144.256419, 144.242852, 144.261894, 144.26052)
wet.dry <- rep(c("dry", "dry", "wet", "wet"),6)

coord <- data.frame(site, site.plot, lat, long, wet.dry)

Summarised the total number of species per method for each plot across all sites.

reptile_summary <- reptile_data %>% 
  unite(site.plot, c(site, plot), sep = ".", remove = FALSE) %>% 
  select(-plot) %>% 
  group_by(site, site.plot, assessment.method) %>% 
  summarise(richness = length(unique(scientific.name))) %>% 
  ungroup() %>% 
  add_column(wet.dry = ifelse(.$site.plot %in% c("Duval.WetA", "Duval.WetB",
                                                 "Tarcutta.WetA", "Tarcutta.WetB",
                                                 "Mourachan.WetA", "Mourachan.WetB",
                                                 "Wambiana.WetA", "Wambiana.WetB",
                                                 "Undara.WetA", "Undara.WetB",
                                                 "Rinyirru.WetA", "Rinyirru.WetB"), "wet", "dry"))

The data frame was missing zero values for methods that didn’t capture any species at a given plot. The following adds zeros for those instances.

# Created object with all combinations of plots and methods to include zero values
plots <- unique(reptile_summary$site.plot)
methods <- unique(reptile_summary$assessment.method)
zeros <- crossing(plots, methods)

# Added reference column for future join 
zeros$plots.methods <- paste(zeros$plots, zeros$methods)

# Selected only columns relevant to join
reptile_summary2 <- reptile_summary %>% 
  select("site.plot", "assessment.method", "richness")

# Added reference column for future join 
reptile_summary2$plots.methods <- paste(reptile_summary2$site.plot, reptile_summary2$assessment.method)

# Joined objects to include zero values for methods across plots
reptile_summary_all <- merge(reptile_summary2, zeros, by = "plots.methods",
                             all.x = T, all.y = T) %>% 
  select("plots", "methods", "richness") %>% 
  replace(is.na(.), 0) %>% 
  rename(site.plot = plots, assessment.method = methods) %>% 
  left_join(coord, by = "site.plot") %>% 
  mutate(site = factor(site, levels = c("Tarcutta", "Duval", "Mourachan", "Wambiana",
                                        "Undara", "Rinyirru")),
         site.plot = factor(site.plot, levels = c("Tarcutta.DryA", "Tarcutta.DryB", "Tarcutta.WetA", "Tarcutta.WetB", "Duval.DryA", "Duval.DryB", "Duval.WetA", "Duval.WetB", "Mourachan.DryA", "Mourachan.DryB", "Mourachan.WetA", "Mourachan.WetB","Wambiana.DryA", "Wambiana.DryB", "Wambiana.WetA", "Wambiana.WetB", "Undara.DryA", "Undara.DryB", "Undara.WetA", "Undara.WetB", "Rinyirru.DryA", "Rinyirru.DryB", "Rinyirru.WetA", "Rinyirru.WetB")),
         assessment.method = factor(assessment.method, levels = c("pitfall", "funnel","incidentals","spotlighting","cover board","camera")),
         wet.dry = factor(wet.dry),
         richness = as.numeric(richness))

4.2 Bayesian Analysis

4.2.1 Priors

Defining priors for Bayesian models.

Priors for the intercept.

reptile_summary_all %>% summarise(median(log(richness)))
##   median(log(richness))
## 1              1.098612
#1.1
reptile_summary_all %>% summarise(mad(log(richness)))
##   mad(log(richness))
## 1           1.454177
#1.5

Priors for the slope.

log(sd(reptile_summary_all$richness))/apply(model.matrix(~assessment.method*lat, data = reptile_summary_all), 2, sd)
##                       (Intercept)           assessment.methodfunnel 
##                               Inf                         3.6955179 
##      assessment.methodincidentals     assessment.methodspotlighting 
##                         3.6955179                         3.6955179 
##      assessment.methodcover board           assessment.methodcamera 
##                         3.6955179                         3.6955179 
##                               lat       assessment.methodfunnel:lat 
##                         0.1921267                         0.1433234 
##  assessment.methodincidentals:lat assessment.methodspotlighting:lat 
##                         0.1433234                         0.1433234 
##  assessment.methodcover board:lat       assessment.methodcamera:lat 
##                         0.1433234                         0.1433234
#3.7

4.2.2 Fit model

4.2.2.1 Model 1

priors1 <- prior(normal(1.1,1.5), class  = "Intercept") +
  prior(normal(0,3.7), class = "b") +
  prior(student_t(3,0,3.7), class = "sd")

methods.form1 <- bf(richness ~ assessment.method*scale(lat, scale = FALSE) + (1|site/site.plot),
                 family = poisson(link = "log"))

methods.brm1 <- brm(methods.form1,
                  data = reptile_summary_all,
                  prior = priors1,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 8,
                  warmup = 1000,
                  backend = 'cmdstanr',
                  save_pars = save_pars(all = TRUE))
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
4.2.2.1.1 Predictions
methods.brm1 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.brm1 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

4.2.2.1.2 Prior vs Posterior
methods.brm1b <- methods.brm1 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 8,
         cores = 3,
         backend = "cmdstanr",
         save_pars = save_pars(all = TRUE))
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 2.1 seconds.
## Chain 3 finished in 2.1 seconds.
## Chain 2 finished in 3.3 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.5 seconds.
## Total execution time: 3.4 seconds.
4.2.2.1.3 Predictions
methods.brm1b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.brm1b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

4.2.2.2 Model 2

priors2 <- prior(normal(1.1,1.5), class  = "Intercept") +
  prior(normal(0,3.7), class = "b") +
  prior(student_t(3,0,3.7), class = "sd")

methods.form2 <- bf(richness ~ assessment.method*scale(lat, scale = FALSE) + (1|site/site.plot),
                 family="negbinomial")

methods.brm2 <- brm(methods.form2,
                  data = reptile_summary_all,
                  prior = priors2,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 8,
                  warmup = 1000,
                  backend = 'cmdstanr',
                  save_pars = save_pars(all = TRUE))
## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 1.9 seconds.
## Chain 3 finished in 3.5 seconds.
## Chain 1 finished in 3.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 3.1 seconds.
## Total execution time: 4.0 seconds.
4.2.2.2.1 Predictions
methods.brm2 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.brm2 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

4.2.2.2.2 Prior vs Posterior
methods.brm2b <- methods.brm2 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 8,
         cores = 3,
         backend = "cmdstanr",
         save_pars = save_pars(all = TRUE))
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 2.9 seconds.
## Chain 3 finished in 2.9 seconds.
## Chain 2 finished in 3.0 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 2.9 seconds.
## Total execution time: 3.1 seconds.
4.2.2.2.3 Predictions
methods.brm2b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.brm2b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

4.2.2.3 Model 3

priors3 <- prior(normal(1.1,1.5), class  = "Intercept") +
  prior(normal(0,3.7), class = "b") +
  prior(student_t(3,0,3.7), class = "sd")

methods.form3 <- bf(richness ~ assessment.method*scale(lat, scale = FALSE) + (1|site/site.plot),
                 family="negbinomial2")

methods.brm3 <- brm(methods.form3,
                  data = reptile_summary_all,
                  prior = priors3,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 1,
                  warmup = 1000,
                  control = list(adapt_delta=0.99),
                  backend = 'cmdstanr',
                  save_pars = save_pars(all = TRUE))
## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 0.5 seconds.
## Chain 1 finished in 0.6 seconds.
## Chain 3 finished in 0.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.7 seconds.
4.2.2.3.1 Predictions
methods.brm3 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.brm3 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

4.2.2.3.2 Prior vs Posterior
methods.brm3b <- methods.brm3 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 2,
         cores = 3,
         control = list(adapt_delta=0.99),
         backend = "cmdstanr",
         save_pars = save_pars(all = TRUE))
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 6.6 seconds.
## Chain 2 finished in 6.6 seconds.
## Chain 3 finished in 9.4 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 7.5 seconds.
## Total execution time: 9.5 seconds.
4.2.2.3.3 Predictions
methods.brm3b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.brm3b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

4.2.2.4 Compare Models

(l.1b <- methods.brm1b %>% loo())
## 
## Computed from 2400 by 144 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -265.3 11.2
## p_loo        16.5  2.2
## looic       530.6 22.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     142   98.6%   427       
##  (0.5, 0.7]   (ok)         2    1.4%   759       
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
(l.2b <- methods.brm2b %>% loo())
## 
## Computed from 2400 by 144 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -266.7 11.1
## p_loo        16.0  2.1
## looic       533.5 22.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     143   99.3%   456       
##  (0.5, 0.7]   (ok)         1    0.7%   660       
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
(l.3b <- methods.brm3b %>% loo())
## 
## Computed from 2400 by 144 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -266.8 11.2
## p_loo        16.2  2.2
## looic       533.5 22.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     141   97.9%   642       
##  (0.5, 0.7]   (ok)         3    2.1%   164       
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
loo_compare(loo(methods.brm1b), loo(methods.brm2b), loo(methods.brm3b))
##               elpd_diff se_diff
## methods.brm1b  0.0       0.0   
## methods.brm2b -1.4       0.4   
## methods.brm3b -1.5       0.4

Model methods.brm1b was selected as best model based on loo estimates.

4.2.3 Diagnostics

methods.brm1b %>% hypothesis("scalelatscaleEQFALSE = 0") %>% plot

4.2.3.1 Trace plots

methods.brm1b$fit %>% stan_trace()

#### Autocorrelation plots

methods.brm1b$fit %>% stan_ac()

#### Rhat statistic

methods.brm1b$fit %>% stan_rhat()

#### Effective sampling size

methods.brm1b$fit %>% stan_ess()

#### Posterior predictive check plot

methods.brm1b %>% pp_check(type = "dens_overlay", ndraws = 200)

#### DHARMa residuals

set.seed(6)
preds <- posterior_predict(methods.brm1b,  ndraws=250,  summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = reptile_summary_all$richness,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = TRUE)
method.resids %>% plot()

#### Dispersion test

method.resids %>% testDispersion()

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.57235, p-value < 2.2e-16
## alternative hypothesis: two.sided

4.2.3.2 Zero-inflation test

method.resids %>% testZeroInflation()

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.1598, p-value = 0.4
## alternative hypothesis: two.sided

4.2.4 Investigation

4.2.4.1 Methods pairwise comparison for each site

newdata_lat <- with(reptile_summary_all,list(lat = c(-35.37876,-30.40221,-27.77876,
                                                     -20.53459,-18.26553,-15.04439)))

(diff.methods3 <- methods.brm1b %>%
  emmeans(~assessment.method|lat, at = newdata_lat) %>%
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n()) %>% #to see if there is any change
  arrange(lat) %>%
  rename("Site" = lat) %>% 
  mutate(Site = as.factor(Site)) %>%  
  mutate(Site = fct_recode(Site, "Tarcutta" = '-35.37876', "Duval" = '-30.40221',
                           "Mourachan" = '-27.77876', "Wambiana" = '-20.53459',
                           "Undara" = '-18.26553', "Rinyirru" = '-15.04439')))
## # A tibble: 90 × 6
## # Groups:   contrast [15]
##    contrast                 Site  Average fractional c…¹ `Lower HDI` `Upper HDI`
##    <fct>                    <fct>                  <dbl>       <dbl>       <dbl>
##  1 pitfall - funnel         Tarc…                   1.05       0.615        1.71
##  2 pitfall - incidentals    Tarc…                   4.04       1.64         7.72
##  3 pitfall - spotlighting   Tarc…                   4.81       2.03         9.65
##  4 pitfall - cover board    Tarc…                   8.03       2.53        18.7 
##  5 pitfall - camera         Tarc…                   8.42       1.79        24.0 
##  6 funnel - incidentals     Tarc…                   3.79       1.59         7.34
##  7 funnel - spotlighting    Tarc…                   4.55       1.67         8.83
##  8 funnel - cover board     Tarc…                   7.58       2.35        17.8 
##  9 funnel - camera          Tarc…                   8.03       1.96        23.6 
## 10 incidentals - spotlight… Tarc…                   1.19       0.356        2.65
## # ℹ 80 more rows
## # ℹ abbreviated name: ¹​`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>

4.2.4.2 Methods pairwise comparison for the average latitude

(diff.meth.avg <- methods.brm1b %>%
  emmeans(~assessment.method|lat) %>% 
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n()) %>% 
  select(-lat))
## # A tibble: 15 × 5
## # Groups:   contrast [15]
##    contrast                   Average fractional chang…¹ `Lower HDI` `Upper HDI`
##    <fct>                                           <dbl>       <dbl>       <dbl>
##  1 pitfall - funnel                                 1.02       0.801        1.28
##  2 pitfall - incidentals                            2.50       1.77         3.37
##  3 pitfall - spotlighting                           3.03       2.12         4.21
##  4 pitfall - cover board                            4.80       2.98         7.09
##  5 pitfall - camera                                11.8        6.74        19.8 
##  6 funnel - incidentals                             2.45       1.79         3.40
##  7 funnel - spotlighting                            2.96       2.04         4.11
##  8 funnel - cover board                             4.67       2.78         6.79
##  9 funnel - camera                                 11.5        5.89        19.0 
## 10 incidentals - spotlighting                       1.21       0.776        1.75
## 11 incidentals - cover board                        1.91       1.13         2.98
## 12 incidentals - camera                             4.69       2.25         7.97
## 13 spotlighting - cover board                       1.58       0.945        2.51
## 14 spotlighting - camera                            3.87       1.87         6.66
## 15 cover board - camera                             2.47       1.08         4.40
## # ℹ abbreviated name: ¹​`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>

4.2.5 Visualisation

Fractional change

methods.em <- methods.brm1b %>%
  emmeans(~assessment.method|lat) %>% 
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(rate = exp(.value))

(methods.mag <- methods.em %>% 
  ggplot() + 
  geom_density_ridges_gradient(aes(x=rate, y=fct_reorder(contrast,rate)),
                               alpha = 0.4, col = "white",
                               quantile_lines = TRUE, quantiles = c(0.025, 0.975),
                               show.legend = FALSE, fill = "#05596E") +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_continuous("Fractional change", trans = "log2",
                     breaks = c(1, 2, 5, 10, 50)) +
  scale_y_discrete(name = "") +
    my.theme())

Proportion of richness across latitudes

lat.list <- with(reptile_summary_all, list(lat = seq(min(lat), max(lat), length = 100)))

newdata <- emmeans(methods.brm1b, ~lat|assessment.method, at=lat.list, type = "response") %>%
  gather_emmeans_draws() %>%
  mutate(richness = exp(.value)) %>%
  group_by(lat, .draw) %>%
  mutate(sum_rate = sum(richness)) %>%
  ungroup() %>%
  mutate(proportion = richness/sum_rate) %>%
  filter(.draw %in% sample(1:2400, 200))

newdata3 <- emmeans(methods.brm1b, ~lat|assessment.method, at=lat.list, type = "response") %>% 
  as.data.frame %>% 
  group_by(lat) %>% 
  mutate(sum_rate.avg = sum(rate),
         proportion.avg = rate/sum_rate.avg) %>% 
  ungroup()

reptile_summary_all_prop <- reptile_summary_all %>% 
  group_by(lat) %>% 
  mutate(sum_rich = sum(richness),
         proportion = richness/sum_rich)


(methods.spag <- newdata %>% 
  ggplot() +
  geom_jitter(data = reptile_summary_all_prop, aes(x = lat, y = proportion, col = assessment.method),
              height = 0, width = 0.2, alpha = 0.4) +
  geom_line(aes(lat, proportion, col = assessment.method, group = interaction(assessment.method,.draw)), alpha=0.07) +
  geom_line(data = newdata3, aes(lat, proportion.avg, group = assessment.method), linewidth = 1.7, col = "black") +
  geom_line(data = newdata3, aes(lat, proportion.avg, col = assessment.method), linewidth = 1) +
  my.theme() +
  scale_y_continuous(name = "Proportion of total species richness",
                     labels = scales::percent,
                     breaks = seq(0, 0.5, by = 0.1),
                     limits = c(-0.025,0.5)) +
  scale_x_continuous(name = "", breaks=c(-35,-30,-25,-20,-15),
                     labels=c("35ºS", "30ºS", "25ºS", "20ºS", "15ºS")) +
  scale_fill_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
                      name = "Sampling Methods", labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
  scale_colour_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
                      name = "", labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
  theme(legend.text = element_text(size = 16, margin = margin(t = 5)),
        legend.position = c(y=0.45, x=-0.07),
        legend.background = element_rect(fill = "transparent")) +
  guides(colour = guide_legend(nrow = 1)) +
  annotate("text", x = c(-35.37876,-30.40221,-27.77876,-20.53459,-18.26553,-15.04439),
           y = c(rep(-0.025, 6)),
           label = c("Tarcutta","Duval","Mourachan","Wambiana","Undara","Rinyirru"), size = 5))

5 Survey Effort

Species accumulation curves

5.1 Prepare data

Merged duplicate rows for duplicate observations and transformed the dataset to a wide format.

reptile_community <- reptile_data %>% 
  select(site, date, assessment.method,scientific.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "scientific.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site)

Add a category for the total richness.

result.i <- vector("list",length(unique(reptile_community$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(reptile_community$site))){ #looped over each study site
  site.i <- reptile_community[reptile_community$site == unique(reptile_community$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$date))) #created an empty list to be filled by each date
  for(a in 1:length(unique(site.i$date))){ #looped over each date
    
    
    date.a <- site.i[site.i$date == unique(site.i$date)[a],] #subset to the date
    
    meta.a <- cbind.data.frame(date.a$site[1],date.a$date[1],c("total"))
    
    date.a <- date.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(date.a))) #calculated column sums
    
    colnames(total.a) <- colnames(reptile_community)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.i <- do.call("rbind.data.frame",result.i)

reptile_community2 <- rbind.data.frame(reptile_community, result.i)

Add a category for the combined pitfall and funnel richness.

result.j <- vector("list",length(unique(reptile_community$site))) #created an empty list that was filled as the loop ran

for(i in 1:length(unique(reptile_community$site))){ #looped over each study site
  site.i <- reptile_community[reptile_community$site == unique(reptile_community$site)[i],] #subset to the site
  
  result.a <- vector("list",length(unique(site.i$date))) #created an empty list to be filled by each date
  for(a in 1:length(unique(site.i$date))){ #looped over each date
    
    
    date.a <- site.i[site.i$date == unique(site.i$date)[a],] %>% 
      filter(assessment.method %in% c("funnel", "pitfall"))
    
    meta.a <- cbind.data.frame(date.a$site[1],date.a$date[1],c("pitfunnel"))
    
    date.a <- date.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
    
    total.a <- cbind.data.frame(meta.a,rbind.data.frame(colSums(date.a))) #calculated column sums
    
    colnames(total.a) <- colnames(reptile_community)
    
    result.a[[a]] <- total.a} #added method results to list
  
  result.j[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

result.j <- do.call("rbind.data.frame",result.j)

reptile_community2 <- rbind.data.frame(reptile_community2, result.j) %>% 
  drop_na()

5.2 Loop

Created a loop that applied specaccum() to all methods across all sites. The output were values for richness of each method at each site over the course of the survey.

# Created a data frame of zeros with 30 rows and the same number of columns as the total number of species.
extra.rows <- data.frame(matrix(NA, nrow = 30, ncol = (ncol(reptile_community2)))) %>% 
  replace(is.na(.), 0) %>% 
  set_names(names(reptile_community2)) %>% 
  select(-(1:3))

result.i <- vector("list",length(unique(reptile_community2$site))) #created an empty list that was filled as the loop ran

set.seed(8) #set seed for reproducibility
for(i in 1:length(unique(reptile_community2$site))){ #looped over each study site
  site.i <- reptile_community2[reptile_community2$site == unique(reptile_community2$site)[i],] #subset to the site
  
  dates.i <- ifelse(unique(site.i$site) %in% c("Tarcutta", "Undara", "Wambiana", "Rinyirru"), 28, 21) #manual addition of dates depending on site to avoid any issues from 0 animals being detected by any method on a given day 
  
  result.a <- vector("list",length(unique(site.i$assessment.method))) #created an empty list to be filled by each method
  for(a in 1:length(unique(site.i$assessment.method))){ #looped over each method
    
    
    method.a <- site.i[site.i$assessment.method == unique(site.i$assessment.method)[a],] #subset to the method
    
    if(nrow(method.a) > 1){ #only proceeded if detections occurred on more than 1 day for specaccum() to work properly
      sampling.a <- method.a$assessment.method[1] #saved the name of the method
      
      method.a <- method.a %>% select(where(is.numeric)) #removed first 3 columns of metadata and only keeps those with numeric data i.e. species columns
      
      if(nrow(method.a)!=dates.i){ #if there were missing days (i.e., days when nothing was found, use the extra.rows object to add rows of zeros)
        method.a <- rbind.data.frame(method.a,extra.rows[1:(dates.i-nrow(method.a)),])}
      
      accum.a <- specaccum(method.a,"random") #calculated data
      
      #created data frame of metadata, richness, and standard deviation
      accum.a <- cbind.data.frame(rep(site.i$site[1],nrow(method.a)),rep(sampling.a,nrow(method.a)),accum.a$richness,accum.a$sd, 1:nrow(method.a))
      colnames(accum.a) <- c("site","assessment.method","richness","sd", "day")
      
      accum.a$lower.ci <- accum.a$richness-qnorm(0.975)*accum.a$sd/sqrt(100) #calculated lower 95% CI
      accum.a$upper.ci <- accum.a$richness+qnorm(0.975)*accum.a$sd/sqrt(100) #calculated upper 95% CI
      
      result.a[[a]] <- accum.a}} #added method results to list
  
  result.i[[i]] <- do.call("rbind.data.frame",result.a) #compressed list into data frame and add to the result.i list
} #end loop

Transformed resulting list to a data frame.

reptiles_wide.result <- do.call("rbind.data.frame",result.i) %>% #compressed result.i list into data frame
mutate(assessment.method = factor(assessment.method, levels = c("pitfall",
                                                                "funnel",
                                                                "incidentals",
                                                                "spotlighting",
                                                                "cover board",
                                                                "camera",
                                                                "total",
                                                                "pitfunnel")),
         site = factor(site, levels = c("Rinyirru", "Undara", "Wambiana", "Mourachan",
                                        "Duval", "Tarcutta")))

5.3 Visualisation

Plot of rarefaction curves for each method split by sites.

(specaccum_reptiles <- ggplot(reptiles_wide.result, aes(x = day, y = richness, col = assessment.method)) + 
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci,fill=after_scale(alpha(colour, 0.3))),
              linetype = 0.1) +
  geom_vline(xintercept=c(7,14,21, 28), linetype="dotted", colour = "black") +
  geom_line(aes(group = assessment.method)) +
  my.theme() +
  facet_wrap(~site) +
  theme(panel.border = element_rect(fill = NA, color = "black", linewidth = 1.5),
        strip.background = element_rect(fill = "lightgrey"),
        axis.title = element_text(size = 26),
        axis.text = element_text(size = 26),
        legend.text = element_text(size = 26),
        legend.position = "bottom",
        strip.text.x = element_text(size = 26)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.07)), breaks = seq(0, 45, by = 5),
                     name = "Species Richness") +
  scale_x_continuous(breaks = seq(0, 28, by = 7),
                     name = "Survey Effort (Days)") +
  scale_colour_manual(values = c("#FF8300","#D14103","#0CC170","black",
                                 "#4E84C4","#8348D8","#FBA6EA","#E8BB5A"),
                      name = "",
                      labels = c("Pitfall Trap","Funnel Trap","Incidental",
                                 "Spotlighting","Arboreal Cover Board", "Camera Trap",
                                 "Total Richness", "Pit + Funnel")))

6 Community composition

6.1 Prepare data

reptile_community3 <- reptile_data %>% 
  filter(!recapture %in% c("y")) %>%  
  select(site, plot, assessment.method,scientific.name) %>% #selected relevant grouping variables
  group_by(across(everything())) %>% #groups by all available columns
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "scientific.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site) %>% 
  mutate(assessment.method = factor(assessment.method))
duv <- reptile_community3 %>% filter(site=="Duval") %>% select(where(~ any(. != 0)))
mou <- reptile_community3 %>% filter(site=="Mourachan") %>% select(where(~ any(. != 0)))
rin <- reptile_community3 %>% filter(site=="Rinyirru") %>% select(where(~ any(. != 0)))
tar <- reptile_community3 %>% filter(site=="Tarcutta") %>% select(where(~ any(. != 0)))
und <- reptile_community3 %>% filter(site=="Undara") %>% select(where(~ any(. != 0)))
wam <- reptile_community3 %>% filter(site=="Wambiana") %>% select(where(~ any(. != 0)))
all <- reptile_community3 %>% select(where(~ any(. != 0)), -"Chelodina longicollis") %>% filter(rowSums(.[,4:ncol(.)]) > 0)
#Function decostand() standardises proportions (method = "total") by rows (MARGIN = 1)
rep_com_bray_duv <- duv %>%
  select(where(is.numeric)) %>%
  decostand(method="total", MARGIN = 1)
rep_com_bray_mou <- mou %>%
  select(where(is.numeric)) %>%
  decostand(method="total", MARGIN = 1)
rep_com_bray_rin <- rin %>%
  select(where(is.numeric)) %>%
  decostand(method="total", MARGIN = 1)
rep_com_bray_tar <- tar %>%
  select(where(is.numeric)) %>%
  decostand(method="total", MARGIN = 1)
rep_com_bray_und <- und %>%
  select(where(is.numeric)) %>%
  decostand(method="total", MARGIN = 1)
rep_com_bray_wam <- wam %>%
  select(where(is.numeric)) %>%
  decostand(method="total", MARGIN = 1)
rep_com_bray_all <- all %>%
  select(where(is.numeric)) %>%
  decostand(method="total", MARGIN = 1)
#Function decostand() standardises presence/absence (method = "pa") by rows (MARGIN = 1)
rep_com_jacc_duv <- duv %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1)
rep_com_jacc_mou <- mou %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1)
rep_com_jacc_rin <- rin %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1)
rep_com_jacc_tar <- tar %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1)
rep_com_jacc_und <- und %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1)
rep_com_jacc_wam <- wam %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1)
rep_com_jacc_all <- all %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1)

6.2 NMDS Ordination

6.2.1 Bray-Curtis dissimilarity

set.seed(1)
all_nmds_bray <-  metaMDS(rep_com_bray_all, distance="bray", k=3,trymax=100, noshare=0.1)
rin_nmds_bray <-  metaMDS(rep_com_bray_rin, distance="bray", k=2,trymax=100)
und_nmds_bray <-  metaMDS(rep_com_bray_und, distance="bray", k=2,trymax=100)
wam_nmds_bray <-  metaMDS(rep_com_bray_wam, distance="bray", k=2,trymax=100)
mou_nmds_bray <-  metaMDS(rep_com_bray_mou, distance="bray", k=2,trymax=1000)

#removed due to insufficient data
duv_nmds_bray <-  metaMDS(rep_com_bray_duv, distance="bray", k=2,trymax=100)
## Warning in metaMDS(rep_com_bray_duv, distance = "bray", k = 2, trymax = 100):
## stress is (nearly) zero: you may have insufficient data
tar_nmds_bray <-  metaMDS(rep_com_bray_tar, distance="bray", k=2,trymax=100)
## Warning in metaMDS(rep_com_bray_tar, distance = "bray", k = 2, trymax = 100):
## stress is (nearly) zero: you may have insufficient data

6.2.2 Jaccard dissimilarity

set.seed(11)
all_nmds_jacc <-  metaMDS(rep_com_jacc_all, distance="jaccard", k=2,trymax=100, noshare=0.1)
rin_nmds_jacc <-  metaMDS(rep_com_jacc_rin, distance="jaccard", k=2,trymax=100)
und_nmds_jacc <-  metaMDS(rep_com_jacc_und, distance="jaccard", k=2,trymax=100)
wam_nmds_jacc <-  metaMDS(rep_com_jacc_wam, distance="jaccard", k=2,trymax=100)
mou_nmds_jacc <-  metaMDS(rep_com_jacc_mou, distance="jaccard", k=2,trymax=100)

#removed due to insufficient data
duv_nmds_jacc <-  metaMDS(rep_com_jacc_duv, distance="jaccard", k=2,trymax=100)
## Warning in metaMDS(rep_com_jacc_duv, distance = "jaccard", k = 2, trymax =
## 100): stress is (nearly) zero: you may have insufficient data
tar_nmds_jacc <-  metaMDS(rep_com_jacc_tar, distance="jaccard", k=2,trymax=100)
## Warning in metaMDS(rep_com_jacc_tar, distance = "jaccard", k = 2, trymax =
## 100): stress is (nearly) zero: you may have insufficient data

6.3 Visualisation

6.3.1 Bray-Curtis dissimilarity

create_plot <- function(data, nmds_data) {
  methods <- nmds_data$points
  species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
  df <- cbind.data.frame(data[,1:3], methods)
  gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
    mutate(assessment.method = factor(assessment.method, levels = c("pitfall","funnel","incidentals","spotlighting","cover board","camera","total")))
  ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
    geom_point(size=3) +
    geom_point(aes(x=mean.x,y=mean.y),size=5) +
    scale_colour_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
                        name = "",
                        labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
    geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) +
    my.theme() + theme(legend.position = "bottom")
}

mds.all.plot.bray.centroid <- create_plot(all, all_nmds_bray)
mds.rin.plot.bray.centroid <- create_plot(rin, rin_nmds_bray)
mds.und.plot.bray.centroid <- create_plot(und, und_nmds_bray)
mds.wam.plot.bray.centroid <- create_plot(wam, wam_nmds_bray)
mds.mou.plot.bray.centroid <- create_plot(mou, mou_nmds_bray)

(mds.bray.all <- (mds.all.plot.bray.centroid + ggtitle('All sites') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.rin.plot.bray.centroid + ggtitle('Rinyirru') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.und.plot.bray.centroid + ggtitle('Undara') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.wam.plot.bray.centroid + ggtitle('Wambiana') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.mou.plot.bray.centroid + ggtitle('Mourachan') & theme(legend.position = "none", plot.title = element_text(size = 24))) + 
   theme(legend.position = c(2,0.65), legend.text = element_text(size = 24)))

6.3.2 Jaccard dissimilarity

create_plot <- function(data, nmds_data) {
  methods <- nmds_data$points
  species <- as.data.frame(nmds_data$species) %>% mutate(species = row.names(.))
  df <- cbind.data.frame(data[,1:3], methods)
  gg <- merge(df, aggregate(cbind(mean.x=MDS1,mean.y=MDS2)~assessment.method, df, mean), by="assessment.method") %>%
    mutate(assessment.method = factor(assessment.method, levels = c("pitfall","funnel","incidentals","spotlighting","cover board","camera","total")))
  ggplot(gg, aes(MDS1,MDS2,color=assessment.method)) +
    geom_point(size=3) +
    geom_point(aes(x=mean.x,y=mean.y),size=5) +
    scale_colour_manual(values = c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8"),
                        name = "",
                        labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
    geom_segment(aes(x=mean.x, y=mean.y, xend=MDS1, yend=MDS2), alpha = 0.2, linewidth = 1.5) +
    my.theme() + theme(legend.position = "bottom")
}

mds.all.plot.jacc.centroid <- create_plot(all, all_nmds_jacc)
mds.rin.plot.jacc.centroid <- create_plot(rin, rin_nmds_jacc)
mds.und.plot.jacc.centroid <- create_plot(und, und_nmds_jacc)
mds.wam.plot.jacc.centroid <- create_plot(wam, wam_nmds_jacc)
mds.mou.plot.jacc.centroid <- create_plot(mou, mou_nmds_jacc)

(mds.jacc.all <- (mds.all.plot.jacc.centroid + ggtitle('All sites') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.rin.plot.jacc.centroid + ggtitle('Rinyirru') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.und.plot.jacc.centroid + ggtitle('Undara') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.wam.plot.jacc.centroid + ggtitle('Wambiana') & theme(legend.position = "none", plot.title = element_text(size = 24))) +
   (mds.mou.plot.jacc.centroid + ggtitle('Mourachan') & theme(legend.position = "none", plot.title = element_text(size = 24))) + 
   theme(legend.position = c(2,0.65), legend.text = element_text(size = 24)))

6.4 Investigation

6.5 Adonis

6.5.1 Bray-Curtis dissimilarity

all_b <- all %>% mutate(assessment.method = factor(assessment.method))
all.dist <- vegdist(all[4:ncol(all)], 'bray')
(all.adonis <- adonis2(all.dist ~ assessment.method, data=all_b))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = all.dist ~ assessment.method, data = all_b)
##                    Df SumOfSqs      R2      F Pr(>F)    
## assessment.method   5   10.309 0.21227 5.5511  0.001 ***
## Residual          103   38.257 0.78773                  
## Total             108   48.566 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Permutation test: significant difference between assessment methods for reptile communities.

#Let's which method is different from which
(adonis.methods.bray <- adonis.pair(all.dist, all$assessment.method))
##                     combination SumsOfSqs   MeanSqs    F.Model         R2
## 1        camera <-> cover board 2.8879066 2.8879066  9.3629077 0.27247135
## 2             camera <-> funnel 2.7328065 2.7328065  7.5459644 0.19081503
## 3        camera <-> incidentals 2.4035118 2.4035118  6.8296694 0.20803345
## 4            camera <-> pitfall 2.7403269 2.7403269  7.4979872 0.18983213
## 5       camera <-> spotlighting 3.1100425 3.1100425 11.2489845 0.31912932
## 6        cover board <-> funnel 2.6722775 2.6722775  6.9302710 0.15088679
## 7   cover board <-> incidentals 1.5522584 1.5522584  4.0656878 0.10968872
## 8       cover board <-> pitfall 2.2871233 2.2871233  5.8898007 0.13120577
## 9  cover board <-> spotlighting 1.0719637 1.0719637  3.2952071 0.09608360
## 10       funnel <-> incidentals 1.5296767 1.5296767  3.7134502 0.08494983
## 11           funnel <-> pitfall 0.3802544 0.3802544  0.9195332 0.01959809
## 12      funnel <-> spotlighting 2.8043239 2.8043239  7.6321637 0.16725404
## 13      incidentals <-> pitfall 1.4649808 1.4649808  3.5336081 0.08116966
## 14 incidentals <-> spotlighting 1.5851769 1.5851769  4.4019111 0.12092527
## 15     pitfall <-> spotlighting 2.7086276 2.7086276  7.3160470 0.16144495
##        P.value P.value.corrected
## 1  0.000999001       0.001152693
## 2  0.000999001       0.001152693
## 3  0.000999001       0.001152693
## 4  0.000999001       0.001152693
## 5  0.000999001       0.001152693
## 6  0.000999001       0.001152693
## 7  0.000999001       0.001152693
## 8  0.000999001       0.001152693
## 9  0.003996004       0.004281433
## 10 0.000999001       0.001152693
## 11 0.502497502       0.502497502
## 12 0.000999001       0.001152693
## 13 0.000999001       0.001152693
## 14 0.000999001       0.001152693
## 15 0.000999001       0.001152693
#comparing all methods with each other
#evidence that everything but pitfall vs funnel traps capture significantly different reptile communities.

6.5.2 Jaccard dissimilarity

all_j <- all %>% mutate(assessment.method = factor(assessment.method))
all.dist2 <- vegdist(rep_com_jacc_all, 'jaccard')
(all.adonis <- adonis2(all.dist2 ~ assessment.method, data=all_j))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = all.dist2 ~ assessment.method, data = all_j)
##                    Df SumOfSqs      R2      F Pr(>F)    
## assessment.method   5    9.069 0.18829 4.7785  0.001 ***
## Residual          103   39.098 0.81171                  
## Total             108   48.168 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Permutation test: significant difference between assessment methods for reptile communities.

#Let's which method is different from which
(adonis.methods.jacc <- adonis.pair(all.dist2, all$assessment.method))
##                     combination SumsOfSqs   MeanSqs    F.Model         R2
## 1        camera <-> cover board 3.0365859 3.0365859 10.3763129 0.29331245
## 2             camera <-> funnel 2.6272329 2.6272329  7.0453213 0.18043958
## 3        camera <-> incidentals 2.3585467 2.3585467  6.5429368 0.20105551
## 4            camera <-> pitfall 2.7649371 2.7649371  7.6180333 0.19228701
## 5       camera <-> spotlighting 2.8379001 2.8379001  9.2684496 0.27859578
## 6        cover board <-> funnel 2.6169850 2.6169850  6.8014219 0.14849805
## 7   cover board <-> incidentals 1.3558707 1.3558707  3.5952112 0.09824267
## 8       cover board <-> pitfall 2.3940005 2.3940005  6.3568896 0.14015268
## 9  cover board <-> spotlighting 1.1861516 1.1861516  3.5283038 0.10218584
## 10       funnel <-> incidentals 1.0036015 1.0036015  2.3528111 0.05555265
## 11           funnel <-> pitfall 0.3502649 0.3502649  0.8346477 0.01782116
## 12      funnel <-> spotlighting 2.0596818 2.0596818  5.2045248 0.12046249
## 13      incidentals <-> pitfall 1.1301393 1.1301393  2.6998874 0.06322938
## 14 incidentals <-> spotlighting 0.5189479 0.5189479  1.3308740 0.03992917
## 15     pitfall <-> spotlighting 2.1228472 2.1228472  5.4802619 0.12604022
##        P.value P.value.corrected
## 1  0.000999001       0.001152693
## 2  0.000999001       0.001152693
## 3  0.000999001       0.001152693
## 4  0.000999001       0.001152693
## 5  0.000999001       0.001152693
## 6  0.000999001       0.001152693
## 7  0.000999001       0.001152693
## 8  0.000999001       0.001152693
## 9  0.000999001       0.001152693
## 10 0.000999001       0.001152693
## 11 0.597402597       0.597402597
## 12 0.000999001       0.001152693
## 13 0.000999001       0.001152693
## 14 0.156843157       0.168046239
## 15 0.000999001       0.001152693
#comparing all methods with each other
#evidence that everything but pitfall vs funnel traps capture significantly different reptile communities.

7 Detection probability

Detection probability between methods for reptile families (maybe also terrestrial vs arboreal)

Binary linear logistic model for the likelihood to detect a reptile of the different families for the methods across sites.

7.1 Prepare data

#Number of individuals of each species captured by each method per survey plot for every survey day.
reptile_community4 <- reptile_data %>% 
  unite(site.plot, c(site, plot), sep = ".", remove = FALSE) %>% 
  select(-plot) %>%
  select(site, site.plot, date, assessment.method, scientific.name) %>% #selected relevant grouping variables
  group_by_at(1:4) %>%
  mutate(number = n()) %>% #created a numbers column that counts the number of duplicate observations
  ungroup() %>% 
  unique() %>% #merged the duplicate rows
  pivot_wider(names_from = "scientific.name", #transformed dataset from long to wide format
              values_from = "number",
              values_fill = list(number=0)) %>% #replaced NAs with 0
  arrange(site)

#Turn abundances for each species into a binary present/absent (1/0).
rep_com_binary <- reptile_community4 %>%
  select(where(is.numeric)) %>%
  decostand(method="pa", MARGIN = 1) %>% 
  cbind(reptile_community4[,c(2,4)]) %>% 
  select(site.plot, assessment.method, everything())

#Make data frame with each method for each site.plot
all.methods <- cbind.data.frame(site.plot=rep(unique(rep_com_binary$site.plot),each=6),asssessment.method=rep(unique(rep_com_binary$assessment.method),24))

#Add merged column
all.methods$site.plot.method <- paste(all.methods$site.plot,all.methods$asssessment.method)

#Find site.plot.methods that are missing from the rep_com_binary because they never documented any species at that plot
missing <- setdiff(all.methods$site.plot.method,paste(rep_com_binary$site.plot,rep_com_binary$assessment.method))

#Make data frame of those missing values with a zero filled in for each species
missing <- cbind.data.frame(all.methods[which(all.methods$site.plot.method %in% missing),1:2],as.data.frame(matrix(0, ncol=ncol(rep_com_binary)-2, nrow=length(missing))))
colnames(missing) <- colnames(rep_com_binary)

#Add missing entries
rep_com_binary <- rep_com_binary %>% 
  bind_rows(missing)
     
#Convert to long format
rep_com_binary.long <- pivot_longer(rep_com_binary, cols = -c(site.plot, assessment.method),
                                    names_to = "scientific.name", values_to = "value")

#Get a count of days when a species was present for a given method for a given plot
days.present <- rep_com_binary.long %>% 
  group_by(site.plot,assessment.method,scientific.name) %>% 
  summarize(days.present = sum(value,na.rm=T))

#Add the survey days information (days surveyed)
survey_days <- read.csv("survey_days.csv")
days.present <- merge(days.present,survey_days,by="site.plot",all.x=T,all.y=T)

#Make a new column that is the total number of days a species was present per plot, summed across all methods and exclude plots where species was never found.
days.present <- days.present %>% 
  group_by(site.plot,scientific.name) %>% 
  mutate(count.per.site.plot = sum(days.present,na.rm=T),
         assessment.method = factor(assessment.method, levels = c("pitfall", "funnel","incidentals","spotlighting","cover board","camera"))) %>% 
  filter(count.per.site.plot > 0) %>% #Remove plots where species were never found
  select(-count.per.site.plot) %>% 
  mutate(det.prob = days.present/days) %>% #Probability of detecting a species via a given method at a given plot.
  ungroup()

This data frame now contains rows for species only for the plots where it was found but for each plot where it was found, it contains rows for each method, regardless whether that method detected it or not.

The days.present column shows the number of days when a species was detected in that plot by that method, the days column shows the number of days when that method was used at that plot and the det.prob column shows the probability of detecting a species via that method at that plot.

Include lifestyle information.

#Join with data frame that has plot and coordinate information
det.plot <- days.present %>% 
  left_join(coord)

#Lifestyles for species per plot
det.plot.groups <- det.plot %>%
  mutate(assessment.method = factor(assessment.method, levels = c("pitfall", "funnel", "incidentals",
                                                                  "spotlighting",
                                                                  "cover board","camera"))) %>%
  add_column(lifestyle = ifelse(.$scientific.name %in% c("Cryptoblepharus australis",
                                                         "Cryptoblepharus pannosus",
                                                   "Diplodactylus vittatus", "Gehyra dubia",
                                                   "Amalosia rhombifer", "Cryptoblepharus adamsi",
                                                   "Cryptoblepharus metallicus", "Oedura castelnaui",
                                                   "Cryptoblepharus virgatus", "Chlamydosaurus kingii",
                                                   "Varanus tristis", "Strophurus williamsi",
                                                   "Oedura coggeri", "Varanus scalaris",
                                                   "Egernia striolata", "Christinus marmoratus",
                                                   "Boiga irregularis", "Hoplocephalus bitorquatus"),
                                "Arboreal Species", "Ground-dwelling Species"))

7.2 Bayesian Analysis - across all methods for all reptiles

7.2.1 Priors

Defining priors for Bayesian models.

Priors for the intercept.

det.plot.groups %>% 
  filter(days.present > 0) %>% 
  summarise(median(days.present))
## # A tibble: 1 × 1
##   `median(days.present)`
##                    <dbl>
## 1                      2
#2
det.plot.groups %>% 
  filter(days.present > 0) %>%
  summarise(mad(days.present))
## # A tibble: 1 × 1
##   `mad(days.present)`
##                 <dbl>
## 1                1.48
#1.5

Priors for the slope.

sd(det.plot.groups$days.present)/apply(model.matrix(~assessment.method*lat, data = det.plot.groups), 2, sd)
##                       (Intercept)           assessment.methodfunnel 
##                               Inf                         7.2485665 
##      assessment.methodincidentals     assessment.methodspotlighting 
##                         7.2485665                         7.2485665 
##      assessment.methodcover board           assessment.methodcamera 
##                         7.2485665                         7.2485665 
##                               lat       assessment.methodfunnel:lat 
##                         0.4340015                         0.3212041 
##  assessment.methodincidentals:lat assessment.methodspotlighting:lat 
##                         0.3212041                         0.3212041 
##  assessment.methodcover board:lat       assessment.methodcamera:lat 
##                         0.3212041                         0.3212041
#4.7 <- too high, let's go for 1.5

7.2.2 Fit model

7.2.2.1 Model 1

priors1 <- prior(normal(2,1.5), class  = "Intercept") +
  prior(normal(0,1.5), class = "b") +
  prior(student_t(3,0,1.5), class = "sd")

methods.form1 <- bf(days.present ~ assessment.method*lat + offset(log(days)) + (1|site/site.plot) +
                      (1|scientific.name), family=poisson(link='log'))

methods.det.brm1 <- brm(methods.form1,
                  data = det.plot.groups,
                  prior = priors1,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 1,
                  warmup = 1000,
                  backend = 'cmdstanr')
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.
7.2.2.1.1 Predictions
methods.det.brm1 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.det.brm1 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

7.2.2.1.2 Prior vs Posterior
methods.det.brm1b <- methods.det.brm1 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 1,
         cores = 3,
         backend = "cmdstanr")
## Running MCMC with 3 parallel chains...
## 
## Chain 3 finished in 128.7 seconds.
## Chain 1 finished in 167.9 seconds.
## Chain 2 finished in 170.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 155.7 seconds.
## Total execution time: 170.6 seconds.
7.2.2.1.3 Predictions
methods.det.brm1b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.det.brm1b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

7.2.2.2 Model 2

priors1 <- prior(normal(2,1.5), class  = "Intercept") +
  prior(normal(0,1.5), class = "b") +
  prior(student_t(3,0,1.5), class = "sd")

methods.form2 <- bf(days.present ~ assessment.method*lat + offset(log(days)) + (1|site/site.plot) +
                      (1|scientific.name), family="negbinomial")

methods.det.brm2 <- brm(methods.form2,
                  data = det.plot.groups,
                  prior = priors1,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 1,
                  warmup = 1000,
                  backend = 'cmdstanr')
## Running MCMC with 3 parallel chains...
## 
## Chain 3 finished in 3.0 seconds.
## Chain 2 finished in 4.3 seconds.
## Chain 1 finished in 4.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 4.0 seconds.
## Total execution time: 4.6 seconds.
7.2.2.2.1 Predictions
methods.det.brm2 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.det.brm2 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

7.2.2.2.2 Prior vs Posterior
methods.det.brm2b <- methods.det.brm2 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 1,
         cores = 3,
         backend = "cmdstanr")
## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 46.3 seconds.
## Chain 1 finished in 46.4 seconds.
## Chain 3 finished in 71.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 54.9 seconds.
## Total execution time: 71.9 seconds.
7.2.2.2.3 Predictions
methods.det.brm2b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.det.brm2b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

7.2.2.3 Model 3

priors1 <- prior(normal(2,1.5), class  = "Intercept") +
  prior(normal(0,1.5), class = "b") +
  prior(student_t(3,0,1.5), class = "sd")

methods.form3 <- bf(days.present ~ assessment.method*lat + offset(log(days)) + (1|site/site.plot) +
                      (1|scientific.name), family="negbinomial2")

methods.det.brm3 <- brm(methods.form3,
                  data = det.plot.groups,
                  prior = priors1,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 1,
                  warmup = 1000,
                  backend = 'cmdstanr')
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.
7.2.2.3.1 Predictions
methods.det.brm3 %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.det.brm3 %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

7.2.2.3.2 Prior vs Posterior
methods.det.brm3b <- methods.det.brm3 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 1,
         cores = 3,
         backend = "cmdstanr")
## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 46.5 seconds.
## Chain 3 finished in 47.4 seconds.
## Chain 1 finished in 48.9 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 47.6 seconds.
## Total execution time: 49.0 seconds.
7.2.2.3.3 Predictions
methods.det.brm3b %>% ggpredict(~assessment.method|lat) %>% plot(add.data = TRUE)

methods.det.brm3b %>% ggpredict(~lat|assessment.method) %>% plot(add.data = TRUE)

7.2.2.4 Compare Models

(m.1b <- methods.det.brm1b %>% loo())
## Warning: Found 13 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## 
## Computed from 2400 by 2070 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -3055.6 132.6
## p_loo       287.6  26.1
## looic      6111.2 265.3
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     2040  98.6%   245       
##  (0.5, 0.7]   (ok)         17   0.8%   78        
##    (0.7, 1]   (bad)         9   0.4%   27        
##    (1, Inf)   (very bad)    4   0.2%   3         
## See help('pareto-k-diagnostic') for details.
(m.2b <- methods.det.brm2b %>% loo())
## Warning: Found 7 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## 
## Computed from 2400 by 2070 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -2150.7  62.1
## p_loo        61.7   5.7
## looic      4301.3 124.2
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     2060  99.5%   379       
##  (0.5, 0.7]   (ok)          3   0.1%   245       
##    (0.7, 1]   (bad)         7   0.3%   56        
##    (1, Inf)   (very bad)    0   0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.
(m.3b <- methods.det.brm3b %>% loo())
## Warning: Found 6 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## 
## Computed from 2400 by 2070 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -2150.5  62.1
## p_loo        61.4   5.5
## looic      4301.0 124.1
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     2058  99.4%   431       
##  (0.5, 0.7]   (ok)          6   0.3%   213       
##    (0.7, 1]   (bad)         6   0.3%   63        
##    (1, Inf)   (very bad)    0   0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.
loo_compare(loo(methods.det.brm1b), loo(methods.det.brm2b), loo(methods.det.brm3b))
## Warning: Found 13 observations with a pareto_k > 0.7 in model
## 'methods.det.brm1b'. It is recommended to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.
## Warning: Found 7 observations with a pareto_k > 0.7 in model
## 'methods.det.brm2b'. It is recommended to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.
## Warning: Found 6 observations with a pareto_k > 0.7 in model
## 'methods.det.brm3b'. It is recommended to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.
##                   elpd_diff se_diff
## methods.det.brm3b    0.0       0.0 
## methods.det.brm2b   -0.2       0.4 
## methods.det.brm1b -905.1      94.1

Model methods.det.brm2b was selected as best model based on loo estimates.

7.2.3 Diagnostics

7.2.3.1 Trace plots

methods.det.brm2b$fit %>% stan_trace()

#### Autocorrelation plots

methods.det.brm2b$fit %>% stan_ac()

#### Rhat statistic

methods.det.brm2b$fit %>% stan_rhat()

#### Effective sampling size

methods.det.brm2b$fit %>% stan_ess()

#### Posterior predictive check plot

methods.det.brm2b %>% pp_check(x="lat", type="intervals")

#### DHARMa residuals

set.seed(5)
preds <- posterior_predict(methods.det.brm2b,  ndraws=250,  summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = det.plot.groups$days.present,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = TRUE)
method.resids %>% plot()

#### Dispersion test

method.resids %>% testDispersion()

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.38588, p-value = 0.048
## alternative hypothesis: two.sided

7.2.3.2 Zero-inflation test

method.resids %>% testZeroInflation()

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99452, p-value = 0.744
## alternative hypothesis: two.sided

7.2.4 Investigation

7.2.4.1 Methods pairwise comparison for the average latitude

(diff.det.meth.avg <- methods.det.brm2b %>%
  emmeans(~assessment.method|lat) %>% 
  regrid(transform = "none") %>% 
  pairs() %>% 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>% 
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n()) %>% 
   select(-lat))
## # A tibble: 15 × 5
## # Groups:   contrast [15]
##    contrast                   Average fractional chang…¹ `Lower HDI` `Upper HDI`
##    <fct>                                           <dbl>       <dbl>       <dbl>
##  1 pitfall - funnel                                1.11        0.832       1.47 
##  2 pitfall - incidentals                           6.36        4.37        8.76 
##  3 pitfall - spotlighting                          3.99        2.69        5.45 
##  4 pitfall - cover board                           5.39        3.68        7.52 
##  5 pitfall - camera                               27.2        15.4        42.3  
##  6 funnel - incidentals                            5.74        3.91        7.89 
##  7 funnel - spotlighting                           3.61        2.48        5.05 
##  8 funnel - cover board                            4.88        3.10        6.71 
##  9 funnel - camera                                24.4        14.2        39.3  
## 10 incidentals - spotlighting                      0.630       0.418       0.889
## 11 incidentals - cover board                       0.848       0.544       1.19 
## 12 incidentals - camera                            4.27        2.36        6.78 
## 13 spotlighting - cover board                      1.35        0.900       1.91 
## 14 spotlighting - camera                           6.81        3.64       10.7  
## 15 cover board - camera                            5.05        2.94        8.04 
## # ℹ abbreviated name: ¹​`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>

7.3 Bayesian Analysis - across all methods for arboreal vs ground-dwelling reptiles

7.3.1 Fit model

7.3.1.1 Model 1

priors1 <- prior(normal(2,1.5), class  = "Intercept") +
  prior(normal(0,1.5), class = "b") +
  prior(student_t(3,0,1.5), class = "sd")

methods.form1 <- bf(days.present ~ assessment.method*lifestyle + offset(log(days)) + (1|site/site.plot) +
                      (1|scientific.name), family="negbinomial")

methods.det.life.brm1 <- brm(methods.form1,
                  data = det.plot.groups,
                  prior = priors1,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 8,
                  warmup = 1000,
                  backend = 'cmdstanr')
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 3.3 seconds.
## Chain 2 finished in 4.0 seconds.
## Chain 3 finished in 3.9 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 3.7 seconds.
## Total execution time: 4.1 seconds.
7.3.1.1.1 Predictions
methods.det.life.brm1 %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)

methods.det.life.brm1 %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)

7.3.1.1.2 Prior vs Posterior
methods.det.life.brm1b <- methods.det.life.brm1 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 8,
         cores = 3,
         backend = "cmdstanr")
## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 35.1 seconds.
## Chain 1 finished in 35.5 seconds.
## Chain 3 finished in 35.8 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 35.4 seconds.
## Total execution time: 35.9 seconds.
7.3.1.1.3 Predictions
methods.det.life.brm1b %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)

methods.det.life.brm1b %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)

7.3.1.2 Model 2

priors2 <- prior(normal(2,1.5), class  = "Intercept") +
  prior(normal(0,1.5), class = "b") +
  prior(student_t(3,0,1.5), class = "sd")

methods.form2 <- bf(days.present ~ assessment.method*lifestyle + offset(log(days)) + (1|site/site.plot) +
                      (1|scientific.name), family="negbinomial2")

methods.det.life.brm2 <- brm(methods.form2,
                  data = det.plot.groups,
                  prior = priors2,
                  sample_prior = "only",
                  refresh = 0,
                  chains = 3, cores = 3,
                  iter = 5000,
                  thin = 5,
                  seed = 8,
                  warmup = 1000,
                  backend = 'cmdstanr')
## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.
7.3.1.2.1 Predictions
methods.det.life.brm2 %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)

methods.det.life.brm2 %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)

7.3.1.2.2 Prior vs Posterior
methods.det.life.brm2b <- methods.det.life.brm2 %>%
  update(sample_prior = "yes",
         refresh = 0,
         seed = 8,
         cores = 3,
         backend = "cmdstanr")
## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 34.4 seconds.
## Chain 3 finished in 36.2 seconds.
## Chain 1 finished in 59.0 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 43.2 seconds.
## Total execution time: 59.0 seconds.
7.3.1.2.3 Predictions
methods.det.life.brm2b %>% ggpredict(~assessment.method|lifestyle) %>% plot(add.data = TRUE)

methods.det.life.brm2b %>% ggpredict(~lifestyle|assessment.method) %>% plot(add.data = TRUE)

7.3.1.3 Compare Models

(m.1b <- methods.det.life.brm1b %>% loo())
## Warning: Found 2 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## 
## Computed from 2400 by 2070 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -1993.7  58.4
## p_loo        66.9   5.8
## looic      3987.4 116.9
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     2060  99.5%   257       
##  (0.5, 0.7]   (ok)          8   0.4%   31        
##    (0.7, 1]   (bad)         1   0.0%   29        
##    (1, Inf)   (very bad)    1   0.0%   13        
## See help('pareto-k-diagnostic') for details.
(m.2b <- methods.det.life.brm2b %>% loo())
## Warning: Found 3 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## 
## Computed from 2400 by 2070 log-likelihood matrix
## 
##          Estimate    SE
## elpd_loo  -1993.1  58.4
## p_loo        65.9   5.5
## looic      3986.3 116.7
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     2062  99.6%   355       
##  (0.5, 0.7]   (ok)          5   0.2%   196       
##    (0.7, 1]   (bad)         2   0.1%   65        
##    (1, Inf)   (very bad)    1   0.0%   36        
## See help('pareto-k-diagnostic') for details.
loo_compare(loo(methods.det.life.brm1b), loo(methods.det.life.brm2b))
## Warning: Found 2 observations with a pareto_k > 0.7 in model
## 'methods.det.life.brm1b'. It is recommended to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.
## Warning: Found 3 observations with a pareto_k > 0.7 in model
## 'methods.det.life.brm2b'. It is recommended to set 'moment_match = TRUE' in
## order to perform moment matching for problematic observations.
##                        elpd_diff se_diff
## methods.det.life.brm2b  0.0       0.0   
## methods.det.life.brm1b -0.5       0.6

Model methods.det.life.brm1b was selected as best model based on loo estimates.

7.3.2 Diagnostics

7.3.2.1 Trace plots

methods.det.life.brm1b$fit %>% stan_trace()

#### Autocorrelation plots

methods.det.life.brm1b$fit %>% stan_ac()

#### Rhat statistic

methods.det.life.brm1b$fit %>% stan_rhat()

#### Effective sampling size

methods.det.life.brm1b$fit %>% stan_ess()

7.3.2.2 DHARMa residuals

set.seed(5)
preds <- posterior_predict(methods.det.life.brm2b,  ndraws=250,  summary=FALSE)
method.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = det.plot.groups$days.present,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = TRUE)
method.resids %>% plot()

#### Dispersion test

method.resids %>% testDispersion()

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.83255, p-value = 0.656
## alternative hypothesis: two.sided

7.3.2.3 Zero-inflation test

method.resids %>% testZeroInflation()

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0084, p-value = 0.648
## alternative hypothesis: two.sided

7.3.3 Investigation

(diff.det.life.meth.rev <- methods.det.life.brm2b %>%
  emmeans(~assessment.method|lifestyle) %>% 
  regrid(transform = "none") %>% 
  pairs() %>% 
  #pairs(reverse = TRUE) %>% to reverse the contrasts 
  gather_emmeans_draws() %>% 
  mutate(f.change = exp(.value)) %>%
  summarise("Average fractional change" = median(f.change),
            "Lower HDI" = HDInterval::hdi(f.change)[1],
            "Upper HDI" = HDInterval::hdi(f.change)[2],
            "Probability of difference" = sum(.value > 0)/n()) %>% 
  arrange(lifestyle))
## # A tibble: 30 × 6
## # Groups:   contrast [15]
##    contrast             lifestyle Average fractional c…¹ `Lower HDI` `Upper HDI`
##    <fct>                <fct>                      <dbl>       <dbl>       <dbl>
##  1 pitfall - funnel     Arboreal…                  3.85       1.70         6.62 
##  2 pitfall - incidenta… Arboreal…                  2.98       1.52         4.95 
##  3 pitfall - spotlight… Arboreal…                  0.541      0.309        0.864
##  4 pitfall - cover boa… Arboreal…                  0.452      0.279        0.718
##  5 pitfall - camera     Arboreal…                 27.9        7.01        71.6  
##  6 funnel - incidentals Arboreal…                  0.770      0.335        1.44 
##  7 funnel - spotlighti… Arboreal…                  0.139      0.0643       0.241
##  8 funnel - cover board Arboreal…                  0.118      0.0551       0.199
##  9 funnel - camera      Arboreal…                  7.26       1.54        20.6  
## 10 incidentals - spotl… Arboreal…                  0.180      0.0984       0.293
## # ℹ 20 more rows
## # ℹ abbreviated name: ¹​`Average fractional change`
## # ℹ 1 more variable: `Probability of difference` <dbl>

7.4 Visualisation

pal <- c("#FF8300", "#D14103","#0CC170","black","#4E84C4","#8348D8")

(det.life.plot <- ggplot(det.plot.groups, aes(x=assessment.method,y=det.prob)) +
  geom_boxplot(aes(color = assessment.method,
                   color = after_scale(darken(color, .1, space = "HLS")),
                   fill = after_scale(desaturate(lighten(color, .8), .4))),
             width = .6, outlier.shape = NA) +
  geom_point(aes(color = assessment.method, color = after_scale(darken(color, .1, space = "HLS"))), fill = "white", shape = 21, stroke = .4, size = 3.5, side = "1",
                            transformation = PositionIdentity) +
  geom_point(aes(fill = assessment.method), color = "transparent", shape = 21, stroke = .4, size = 3.5, alpha = .3, side = "1", transformation = PositionIdentity) +
  facet_wrap(~lifestyle) +
  scale_y_continuous(labels = scales::percent,
                     breaks = seq(0, 0.9, by = 0.1),
                     limits = c(0,0.9),
                     name = "Detection Probability") +
  scale_x_discrete(name = "", guide = "none") +
  scale_fill_manual(values = pal, guide = "none") +
  scale_colour_manual(values = pal, name = "", labels = c("Pitfall Trap","Funnel Trap","Incidental","Spotlighting","Arboreal Cover Board","Camera Trap")) +
  my.theme() +
  theme(panel.grid.major.y = element_line(colour = "grey", linewidth = 0.2, linetype = "dotted"),
        panel.grid.minor.y = element_line(colour = "grey", linewidth = 0.1, linetype = "dotted"),
        legend.position = c(0.5, -0.07),
        legend.background = element_rect(fill = "transparent"),
        legend.direction = "horizontal",
        plot.margin = unit(c(5, 10, 12, 10), units = "mm"),
        panel.border = element_rect(fill = NA, color = "black"),
        strip.background = element_rect(fill = "lightgrey")))
## Warning: Duplicated aesthetics after name standardisation: colour
## Duplicated aesthetics after name standardisation: colour
## Warning in geom_point(aes(color = assessment.method, color =
## after_scale(darken(color, : Ignoring unknown parameters: `side` and
## `transformation`
## Warning in geom_point(aes(fill = assessment.method), color = "transparent", :
## Ignoring unknown parameters: `side` and `transformation`

8 Session Info

sessionInfo()
## R version 4.3.0 (2023-04-21)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.3.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Australia/Brisbane
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] report_0.5.7         HDInterval_0.2.4     patchwork_1.1.2     
##  [4] EcolUtils_0.1        vegan_2.6-4          lattice_0.21-8      
##  [7] permute_0.9-7        ggridges_0.5.4       tidybayes_3.0.4     
## [10] emmeans_1.8.6        DHARMa_0.4.6         ggeffects_1.2.2     
## [13] rstan_2.21.8         StanHeaders_2.21.0-7 brms_2.19.0         
## [16] Rcpp_1.0.10          cmdstanr_0.5.3       colorspace_2.1-0    
## [19] lubridate_1.9.2      forcats_1.0.0        stringr_1.5.0       
## [22] dplyr_1.1.2          purrr_1.0.1          readr_2.1.4         
## [25] tidyr_1.3.0          tibble_3.2.1         ggplot2_3.4.2       
## [28] tidyverse_2.0.0      knitr_1.42          
## 
## loaded via a namespace (and not attached):
##   [1] tensorA_0.36.2       rstudioapi_0.14      jsonlite_1.8.4      
##   [4] magrittr_2.0.3       estimability_1.4.1   farver_2.1.1        
##   [7] nloptr_2.0.3         rmarkdown_2.21       vctrs_0.6.2         
##  [10] minqa_1.2.5          base64enc_0.1-3      htmltools_0.5.5     
##  [13] haven_2.5.2          distributional_0.3.2 sass_0.4.6          
##  [16] bslib_0.4.2          htmlwidgets_1.6.2    plyr_1.8.8          
##  [19] zoo_1.8-12           cachem_1.0.8         igraph_1.4.2        
##  [22] iterators_1.0.14     mime_0.12            lifecycle_1.0.3     
##  [25] pkgconfig_2.0.3      gap_1.5-1            colourpicker_1.2.0  
##  [28] Matrix_1.5-4         R6_2.5.1             fastmap_1.1.1       
##  [31] shiny_1.7.4          digest_0.6.31        ps_1.7.5            
##  [34] qgam_1.3.4           crosstalk_1.2.0      labeling_0.4.2      
##  [37] fansi_1.0.4          timechange_0.2.0     abind_1.4-5         
##  [40] mgcv_1.8-42          compiler_4.3.0       doParallel_1.0.17   
##  [43] withr_2.5.0          backports_1.4.1      inline_0.3.19       
##  [46] shinystan_2.6.0      highr_0.10           pkgbuild_1.4.0      
##  [49] MASS_7.3-58.4        gtools_3.9.4         loo_2.6.0           
##  [52] tools_4.3.0          httpuv_1.6.11        threejs_0.3.3       
##  [55] glue_1.6.2           callr_3.7.3          nlme_3.1-162        
##  [58] promises_1.2.0.1     grid_4.3.0           checkmate_2.2.0     
##  [61] cluster_2.1.4        reshape2_1.4.4       generics_0.1.3      
##  [64] gtable_0.3.3         tzdb_0.4.0           data.table_1.14.8   
##  [67] hms_1.1.3            utf8_1.2.3           foreach_1.5.2       
##  [70] pillar_1.9.0         ggdist_3.3.0         markdown_1.7        
##  [73] posterior_1.4.1      later_1.3.1          splines_4.3.0       
##  [76] tidyselect_1.2.0     miniUI_0.1.1.1       arrayhelpers_1.1-0  
##  [79] gridExtra_2.3        stats4_4.3.0         xfun_0.39           
##  [82] bridgesampling_1.1-2 matrixStats_0.63.0   DT_0.27             
##  [85] stringi_1.7.12       yaml_2.3.7           boot_1.3-28.1       
##  [88] evaluate_0.21        codetools_0.2-19     cli_3.6.1           
##  [91] RcppParallel_5.1.7   shinythemes_1.2.0    xtable_1.8-4        
##  [94] munsell_0.5.0        processx_3.8.1       jquerylib_0.1.4     
##  [97] coda_0.19-4          svUnit_1.0.6         parallel_4.3.0      
## [100] rstantools_2.3.1     ellipsis_0.3.2       prettyunits_1.1.1   
## [103] gap.datasets_0.0.5   dygraphs_1.1.1.6     bayesplot_1.10.0    
## [106] Brobdingnag_1.2-9    lme4_1.1-33          mvtnorm_1.1-3       
## [109] scales_1.2.1         xts_0.13.1           insight_0.19.1      
## [112] crayon_1.5.2         rlang_1.1.1          shinyjs_2.1.0

9 Cite Packages

cite_packages()
##   - Bürkner P (2017). "brms: An R Package for Bayesian Multilevel Models Using Stan." _Journal of Statistical Software_, *80*(1), 1-28. doi:10.18637/jss.v080.i01 <https://doi.org/10.18637/jss.v080.i01>. Bürkner P (2018). "Advanced Bayesian Multilevel Modeling with the R Package brms." _The R Journal_, *10*(1), 395-411. doi:10.32614/RJ-2018-017 <https://doi.org/10.32614/RJ-2018-017>. Bürkner P (2021). "Bayesian Item Response Modeling in R with brms and Stan." _Journal of Statistical Software_, *100*(5), 1-54. doi:10.18637/jss.v100.i05 <https://doi.org/10.18637/jss.v100.i05>.
##   - Eddelbuettel D, François R (2011). "Rcpp: Seamless R and C++ Integration." _Journal of Statistical Software_, *40*(8), 1-18. doi:10.18637/jss.v040.i08 <https://doi.org/10.18637/jss.v040.i08>. Eddelbuettel D (2013). _Seamless R and C++ Integration with Rcpp_. Springer, New York. doi:10.1007/978-1-4614-6868-4 <https://doi.org/10.1007/978-1-4614-6868-4>, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta JJ (2018). "Extending extitR with extitC++: A Brief Introduction to extitRcpp." _The American Statistician_, *72*(1), 28-36. doi:10.1080/00031305.2017.1375990 <https://doi.org/10.1080/00031305.2017.1375990>.
##   - Gabry J, Češnovar R (2022). _cmdstanr: R Interface to 'CmdStan'_. https://mc-stan.org/cmdstanr/, https://discourse.mc-stan.org.
##   - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate." _Journal of Statistical Software_, *40*(3), 1-25. <https://www.jstatsoft.org/v40/i03/>.
##   - Hartig F (2022). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. R package version 0.4.6, <https://CRAN.R-project.org/package=DHARMa>.
##   - Kay M (2023). _tidybayes: Tidy Data and Geoms for Bayesian Models_. doi:10.5281/zenodo.1308151 <https://doi.org/10.5281/zenodo.1308151>, R package version 3.0.4, <http://mjskay.github.io/tidybayes/>.
##   - Lenth R (2023). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. R package version 1.8.6, <https://CRAN.R-project.org/package=emmeans>.
##   - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." _Journal of Open Source Software_, *3*(26), 772. doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
##   - Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023). "Automated Results Reporting as a Practical Tool to Improve Reproducibility and Methodological Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
##   - Meredith M, Kruschke J (2022). _HDInterval: Highest (Posterior) Density Intervals_. R package version 0.2.4, <https://CRAN.R-project.org/package=HDInterval>.
##   - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version 3.2.1, <https://CRAN.R-project.org/package=tibble>.
##   - Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, O'Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J (2022). _vegan: Community Ecology Package_. R package version 2.6-4, <https://CRAN.R-project.org/package=vegan>.
##   - Pedersen T (2022). _patchwork: The Composer of Plots_. R package version 1.1.2, <https://CRAN.R-project.org/package=patchwork>.
##   - R Core Team (2023). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>.
##   - Salazar G (2023). _EcolUtils: Utilities for community ecology analysis_. R package version 0.1, <https://github.com/GuillemSalazar/EcolUtils>.
##   - Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer, New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
##   - Simpson G (2022). _permute: Functions for Generating Restricted Permutations of Data_. R package version 0.9-7, <https://CRAN.R-project.org/package=permute>.
##   - Stan Development Team (2020). "StanHeaders: Headers for the R interface to Stan." R package version 2.21.0-6, <https://mc-stan.org/>.
##   - Stan Development Team (2023). "RStan: the R interface to Stan." R package version 2.21.8, <https://mc-stan.org/>.
##   - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag New York. ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>.
##   - Wickham H (2022). _stringr: Simple, Consistent Wrappers for Common String Operations_. R package version 1.5.0, <https://CRAN.R-project.org/package=stringr>.
##   - Wickham H (2023). _forcats: Tools for Working with Categorical Variables (Factors)_. R package version 1.0.0, <https://CRAN.R-project.org/package=forcats>.
##   - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
##   - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar of Data Manipulation_. R package version 1.1.2, <https://CRAN.R-project.org/package=dplyr>.
##   - Wickham H, Henry L (2023). _purrr: Functional Programming Tools_. R package version 1.0.1, <https://CRAN.R-project.org/package=purrr>.
##   - Wickham H, Hester J, Bryan J (2023). _readr: Read Rectangular Text Data_. R package version 2.1.4, <https://CRAN.R-project.org/package=readr>.
##   - Wickham H, Vaughan D, Girlich M (2023). _tidyr: Tidy Messy Data_. R package version 1.3.0, <https://CRAN.R-project.org/package=tidyr>.
##   - Wilke C (2022). _ggridges: Ridgeline Plots in 'ggplot2'_. R package version 0.5.4, <https://CRAN.R-project.org/package=ggridges>.
##   - Xie Y (2023). _knitr: A General-Purpose Package for Dynamic Report Generation in R_. R package version 1.42, <https://yihui.org/knitr/>. Xie Y (2015). _Dynamic Documents with R and knitr_, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, <https://yihui.org/knitr/>. Xie Y (2014). "knitr: A Comprehensive Tool for Reproducible Research in R." In Stodden V, Leisch F, Peng RD (eds.), _Implementing Reproducible Computational Research_. Chapman and Hall/CRC. ISBN 978-1466561595.
##   - Zeileis A, Fisher JC, Hornik K, Ihaka R, McWhite CD, Murrell P, Stauffer R, Wilke CO (2020). "colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes." _Journal of Statistical Software_, *96*(1), 1-49. doi:10.18637/jss.v096.i01 <https://doi.org/10.18637/jss.v096.i01>. Zeileis A, Hornik K, Murrell P (2009). "Escaping RGBland: Selecting Colors for Statistical Graphics." _Computational Statistics & Data Analysis_, *53*(9), 3259-3270. doi:10.1016/j.csda.2008.11.033 <https://doi.org/10.1016/j.csda.2008.11.033>. Stauffer R, Mayr GJ, Dabernig M, Zeileis A (2009). "Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations." _Bulletin of the American Meteorological Society_, *96*(2), 203-216. doi:10.1175/BAMS-D-13-00155.1 <https://doi.org/10.1175/BAMS-D-13-00155.1>.